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ABSTRACT 


As part of the existing acoustic transient localization program, a feasibility study 
was performed to apply existing algorithms to signals at higher carrier frequencies. 
The coherent matching, autocorrelation matching and SIFT algorithms are time- 
domain Matched Field Processing algorithms based on arrival structures for single 
hydrophone applications. In previous studies, these algorithms were employed only 
at lower frequencies using ray propagation models to create the replicas with 
varying success. This study is meant to investigate the performance of the 
algorithms at higher frequencies, using both the University of Miami Parabolic 
Equation (UMPE) Model and the Hamiltonian Raytracing Program for the Ocean 
(HARPO), to give insight into the previously unexplained inconsistent behavior of 
the algorithms at low frequencies, to improve and optimize existing algorithms, to 
point out improvements to existing eigenray extraction programs, and to suggest 
additional signal processing on the signal. Simulations are performed and synthetic 
signals are generated using both the HARPO and UMPE models. The arrival 
structures are investigated and the relation between features in the arrival structures 
for matching and the physical parameters are identified. Some insight into the 
performance of the SIFT algorithm is gained which relates matching and physical 
parameters. Simulations lead to improvements and optimization of the algorithms 
and give insight into the performance at higher frequencies. 
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I. INTRODUCTION 


Acoustic source localization has been an intensive research issue for the past few 
decades. Prior to that, passive SONAR was used to obtain an estimate of the source direction 
and other techniques, such as Target Motion Analysis (TMA), Contact Motion Analysis 
(CMA) and triangulation, were used to get an estimate of the range to the source. With the 
introduction of faster computers, other possibilities arose. One of the techniques that was 
developed is known as Matched Field Processing (MFP). The MFP process consists of 
systematically placing a synthetic source at each point in a search grid and comparing the 
signal received from all the synthetic sources with the signal received from the true source. 
When the synthetic source location and the tme source location match, the correlation of the 
true and the synthetic received signals should be maximum. Most of the work in this field 
has been done with array receivers and narrow band signals. The project-in-hand concerns 
itself with the case of a single receiver hydrophone and broadband signals. 

The objective of this study is an examination of the influence of the physics 
mismatch in the prediction of the acoustic propagation on a number of MFP algorithms for 
a single hydrophone receiver and a transient-like point source. A time-domain signal 
autocorrelation matching is performed to produce the ambiguity surface for the localization. 
A full wave, parabolic equation model is employed to produce a synthetic signal and a 
reciprocal prediction to provide a baseline for the correlation. Predictions from a ray-based 
propagation model are then matched to the synthetic signal. By comparing this ambiguity 


surface with the baseline, the influence of the physics mismatch due to the ray predictions 


can be quantified. Several aspects of the ray model can be independently affected providing 
information of model degradation on localization performance. A secondary but not less 
important objective is to suggest directions for future research beyond the present available 
algorithms. 

The remainder of this thesis consists of five chapters. Chapter II describes the 
propagation models and arrival structure synthesis. Chapter III describes the matching 
algorithms and the influence of physical parameters on these algorithms. Chapter IV 
describes the setup of the experiments and the results of the autocorrelation matching and 
the approximate autocorrelation matching. Chapter V describes the issues which were 
encountered during the research and were not investigated, but which may lead to more 
insight into the performance of present algorithms or lead to more robust algorithms. 
Chapter VI presents the conclusions of the study. In the Appendix, a derivation for a 
frequency based analogue of a generalized beamformer for a single hydrophone is presented, 


and the relationship to autocorrelation matching 1s briefly discussed. 


Il. PROPAGATION MODELS AND ARRIVAL STRUCTURE SYNTHESIS 


Propagation models and in particular the relation of their parameters to physical 
parameters form a major part of the basis of this research. The step from the output of the 
propagation models to an arrival structure from another important part. This is also the step 
where some of the major approximations are made. For this research, two propagation 
models have been used: a full wave parabolic equation model and a Hamiltonian raytracing 
propagation model. The models are based on the same wave equation but different 
assumptions are made to generate an approximate solution. We begin with the homogeneous 
wave equation for the acoustic pressure in a medium with sound speed c(x) and density 
o(x), (Jensen et al. 1994) 


Vp(x,t)|__ 1 a op%,t) t) 
vy) faces = 
i | p(x) |. esa er? ae () 


Assuming a harmonic solution, p(x,f)=p(x)e"’, leads to the corresponding Helmholtz 


equation 


acar| a ee) OPO) «5.6 (2) 


px) ] c(X) 


where S is the source function and p is pressure, which are generally complex. The nght- 


hand term is zero when the source is not included, which is true for the major part of the 


ocean. For a fixed density and a range independent environment the formulas could even be 


more simplified to 





Sere | 0-0 (3) 
(7) 


Although the models incorporate more features and are based on generalized versions of the 


above formulas (e.g., including currents) those features are not used in this project. 
A. HAMILTONIAN RAYTRACING 


Raytracing or geometrical acoustics assumes that sound waves travel along 
geometric paths called rays. Each ‘ray’ travels along a certain path which is dictated by the 
local sound speed and, where appropriate, boundary interactions. Assuming a constant 
density, we assume a solution of the form (Jensen et al., 1994) 


iWT(x) 
on aa ( - 4) 


which represents a ray series solution. This form, although in general divergent, provides an 
asymptotic approximation to the exact solution of the Helmholtz equation in certain cases. 


The series term expresses the frequency dependence of the amplitude. 


Usually only the first term in the series solution is used in the derivation. The 


goveming equations for the ray theory are then 





ve(xP=—_, 
T(x) Ww (5) 
2V1(X):VA,(X)+V't(X)A=0, (6) 
and 
2WV1(£)'VA (X)tV'O(X)A,=-WA(Z)  f=1,2......, (7) 


where Tis the travel time phase and 4 is the amplitude. The first equation is known as the 
eikonal equation. It determines the local direction of the phase front and, as rays travel 
perpendicular to phase fronts, it determines the ray paths. The second and latter equations 
are called the transport equations of which only the first is commonly used. This is the same 
as taking a first order approximation instead of a series solution. The second and latter 
equations determine the amplitude along the ray path. Furthermore, the expressions above 


make use of the approximations 


2 
w»Vc(x) and =e} 


0 


such that variations of the sound speed and variations in amplitude occur over much larger 
scales than a wavelength. The neglect of higher order terms in the above solution leads to 
the assumption that the energy is conserved within a ray tube (higher order terms include 
leakage) and leads to a break down of the transport equation at or near focal points 
(caustics). 

In geometric acoustics, a ray behaves much like a particle. Therefore, Hamiltonian 
theory can be applied to the propagation of sound (Lighthill, 1978). The first order ray 


equations derived from the eikonal equation can be written as 


ax LS 
ane? K(X), (9) 
dk 
wg 
— ' (11) 


where s 1s the path length, Oe 1s the wave number, i 1s the reference wave 
C(x C 
0 


number, cy, 1s the reference sound speed and n is the acoustical index of refraction. The 


equivalent Hamiltonian equations can be written as 


ax a 
ra ZH(x,k,w,), (12) 


dk 


—=-VH(x,kw,0), 3 
oF ( ) (13) 
H(¥,k,w,t)=w* -c?(*)k?=0, (14) 
and 
al 15 
Anf (15) 


In this information, the Hamiltonian is zero along the ray path. For this project, time and 
frequency dependence of the Hamiltonian are neglected (i.e., no volume attenuation, 
currents, or changing environments). All frequency dependent features such as boundary 
interactions are calculated in the postprocessing. Although the Hamiltonian 1s a constant of 
the motion in this formulation, the number of degrees of freedom (3) exceeds the number 
of constants. In general, ray paths are then chaotic (Smith et al., 1992). However, in this 


project, the environment is two dimensional and range-independent, and all the ray paths are 


regular. 


B. HAMILTONIAN RAYTRACING PROGRAM FOR THE OCEAN (HARPO) 


The raytracing program used for this thesis is the Hamiltonian Raytracing Program 


for the Ocean (HARPO). HARPO is written in Fortran and its roots go back to a program 


written by Dudrak (1961). It has since been updated by many different people. The program 
treats the ocean as a continuous changing environment and not as a large number of 
stratified segments thereby avoiding effects such as false caustics. Of the ray parameters, 
only the ray paths and the travel times are computed by the HARPO program. Amplitude and 
path length are determined in the postprocessing. 

The heart of the program is a fourth-fifth order Adams Moulton predictor corrector 
algorithm to integrate the Hamiltonian equations along the ray path. The accuracy of the 
integration process and the spacing of the initial ray fan are extremely important parameters 
for the calculation of the eigenray parameters. Comparisons with analytical results of travel 
times have indicated that the calculations are highly accurate with phase errors < 1/100 up 
to 10 km at 1 kHz. 

Two features in the design of the program make the standard output not particularly 
suitable for a straightforward eigenray extraction routine. The rays are referenced to a 
spherical coordinate system and the points along the ray are not equally spaced in range. 
This hampers some of the necessary improvements for eigenray extractions close to the 
boundaries, and makes it necessary to use a local Cartesian coordinate system to avoid 
complicated formulations. In what follows, a two-dimensional environment, specifying 


range and depth positions (r,z) will be assumed. 


oF EIGENRAY EXTRACTION 


Eigenrays are rays that connect source and receiver. The number of eigenrays is 
generally quite large. Usually only a small number has enough amplitude to be important, 
however. Eigenray extraction programs are, in a general sense, root solvers. In Fig. 1, a plot 
of launch angle versus depth at some finite range, it is clearly shown that the eigenrays for 


a given range and depth are the roots for that depth of the curve given in the figure. 


depth (m) 


60 
50 
40 
30 





launch angle 


Figure 1 Launch angle versus depth plot at 
a certain range. 


As the output of HARPO is not organized in a convenient range/depth grid, the process of 
finding the eigenrays is more involved. An eigenray extraction program developed in Matlab 
(ray3d.m) by Chiu deals with this matter. The program selects the rays with a (local) 


minimum vertical distance to the receiver location out of the total fan of rays and calculates 


the ray parameters by a third order polynomial interpolation using these selected rays and 
adjacent rays. No retracing of the eigenrays is done. The eigenray parameters calculated by 
the program are: phase travel time 7, the path length s, the local sound speed c, the local 
sound speed gradient %, the local direction of the ray 6, the launch angle @,, the amplitude 
due to reflections, the phase shift due to reflections and turning points, and the amplitude 
due to geometrical spreading. 

The amplitude along the ray is not calculated in HARPO, but in the eigenray 
extraction routine based on the equivalent WKB formulation of the ray solution. The 
calculation is based on the conservation of power within a ray tube. The formula based on 
the assumption of a point source, used for the amplitude calculation is (Ziomek, 1995) 


p,(r.z)c(7,2) R,cos(B,) 


oe 16 
Pol712,)e(7,.2,) risin(B(r,2))dr1/aB a 


P{rz)P=P Ar z)P 


where |P{7,z)| 1s the pressure amplitude at the gridpoint, PAT, Z,)| is the pressure amplitude 
at 1 m (which can be related to the source level), ,/r,z) is the density at the gridpoint, 
P(r ,Z,) 1s the density at 1 m from the source, c(r,z) is the sound speed at the gridpoint, 
c(r,z, is the reference sound speed at 1 m, R,=1 m, f, is the launch angle of the eigenray, 
r is the horizontal range from the source of the gridpoint, and dr relates to the distance 


perpendicular to the direction of the wavefront. We can rewrite this using 


10 


ae hat 
PArzP 
Pea) =Pg\7.Z), (17) 
Pee 


6=rsin(B), 


and arrive at 


e(r,z) ©os(By) 
c(r,z,) rd/dB, 


Q= 








(18) 


This formulation is known to produce infinite amplitudes at and near turning and 
focal points. However, the numerical implementation used here prevents the amplitude from 
blowing up. Ordinarily the 6/d, is calculated using the selected ray (which is locally closest 
to the gridpoint). When this ray is too close, the fan is opened up a little. The same is done 
when the amplitude exceeds the cylindrical spreading amplitude. 

The routines for the reflection coefficients are based on plane wave reflection 
coefficients. The surface reflection coefficient 1s set to unity, with a 180 degrees phase shift, 
R,=e '", while the bottom reflection coefficient is calculated including the sediment 


attenuation « (dB/Hz/km) by 


11 


Z 





cf 
ae 20log,,(e), (19) 
“<p F000 BO? 
ee 
0 | (1) 050) (20) 
Cc 


SO 


_ Pp, (Ce, tic ,,)sin(9,)-pesin(9,) 





ci p,(c, ric, ,sin(O;)+pcsin(O an ea 
Written in terms of an amplitude and a phase this leads to 
R,=|R,le'**, (22) 
where 
c,sin(0.)-pcsin(8,))*+(p,c, sin(8,))’ 
RJ- (P,c, 3 (O)=pesin(,)Y (Pse,,;8iNO)) (23) 
(p,c,sin(8,)+pcsin(6,))’ +(p,¢, sin(®,))’ 
and 
: c, sin(6. sin(8. 
jpctan' P2810) gs PSO) 
p,c,sin(8.)-pcsin(O,) p,c,sin(0,)+pcsin(8,) 
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At ray turning points, we can express the phase change in terms of an equivalent “reflection 
coefficient”, 1.¢.,R, =e "2 The total phase for selected rays (which are representative of 


the eigenrays) is then calculated as 


N, 
0=-N.w2-N.n-)) O, » (25) 
jf : 


where the number of turning points, surface reflections and bottom reflections are given by 
N,, N, and N, respectively. 

The phase of a ray changes according to reflections, refractions and travel time. 
Refractions close to the surface and close to the bottom may not cause a discrete -7t/2 phase 
change but something between -1/2 and the phase change due to the reflection. Furthermore, 
the position in range is not strictly localized at the refracted points. This together with the 
fact that small errors accumulate per bounce and may give the phase of a ray a random 
component which can be quite large. 

The estimation of the amplitude of a ray 1s also numerically difficult. Opening up the 
ray fan as discussed before does not give you a true estimate but rather a conservative 
estimate due to the averaging caused by opening up the ray fan. However, travel times are 
estimated very accurately. Still, at high frequencies, small errors may lead to unacceptable 
phase errors. For ranges and frequencies considered 1n this project, we do not expect travel 


time errors to be significant. 


D. ARRIVAL STRUCTURE SYNTHESIS FROM EIGENRAY DATA 


Several different approaches have been used to derive arrival structures from 
eigenray data for different localization algorithms. The important parameters used in the 


synthesis are 


rs eigenray travel time (relative and absolute travel time), 

0 eigenray phase (due to reflections and refractions), 

fe carner frequency (which together with the travel time may account 
for a major part of the phase), 

as eigenray amplitude, and 

s(t),S(f) source function (pulse/transient) in time or frequency domain (real 


or complex). 
For the synthesis, the ocean is considered as a linear time invariant filter causing a time 
delay and a phase shift. The total amplitude and phase contribution due to the reflections and 


refractions for each eigenray can be written as 


RQH= 








N b N s N t 
I] RM) I] RW) I] RM) 
i p= i 


=IR (fle, 


(26) 


where the product in Eq. (25) represents the accumulated effects of multiple reflections and 


refractions and the last equation specifies the conjugation for negative frequencies. For this 
a+ 6) feo 27) 


project we assume that R, 1s frequency independent. The frequency response of a real 


received signal from an individual arrival can then be written as 


Y Q=7,D) XY (28) 


where 


HA)=4,7) R, 7,2) (29) 


is the complex frequency domain transfer function, a, represents the geometrical amplitude 
factor, and 7, represents the phase change due to the time delay. Note that Y (/)=Y, (-/) 
since the time domain signal is real. We assume that the amplitude factor ts also frequency 


independent and the delay factor can be written as 


G@ise (aay (30) 


15 


We can now write Y,(f) as 


Y (fza, Re me PM K(f), (31) 
In the time domain, this can be written as 
y,(t) =a, R, [Me ts a oP (32) 


The arrival structure is then defined as the sum over all eigenrays of the individual arnvals 


Lom 


N 
WO=d, y,(t). (33) 


For some applications, it may be easier to first sum in the frequency domain, 


Y=) YA). | (34) 


and then transform the results to the time domain, 


Wd= | Ypeaf (35) 
These two approaches should be equivalent. 
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For some localization algorithms, it may be preferable to perform all the calculations 
in the frequency domain. For bandpass signals, it might be desirable to cast it in terms of the 
preenvelope or complex envelope of the received signal. This results in a few minor changes 
in the formulas. Specificaily, we treat the time domain signal as complex (1.e., neglect the 


negative portion of the spectrum) and define the basebanded frequency response as 
VL O=H, Ff) XD) (36) 
or 
Lee a X(f) e ire | Memes = | (37) 
If we write the complex source function as 


RD=F Xf} (38) 


where 


RN =fx(1)+iX(D}e 7" (39) 


then the complex received signal of a single eigenray can be wmitten as 
mp =(O0-=/2Kf08 
WO)=x(t-T Je 


(40) 
=X(t-T )cos(O+27f t_)-ix(t-t)sin(O+277 Tt, ) 


where /, is the carrier frequency or center frequency of the pass band. 
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Another way to write the received signal, which was used in previous simulations 
and which is especially applicable to low pass signals, can be derived starting from the 


preenvelope (no base banding) of the source function 


X(t) =x(t) + x(2), (41) 


where <(f) is the Hilbert transform of the real source function x(t), and the frequency 
spectrum of the preenvelopex (7) is the same as that of x(t) but with the negative frequencies 


discarded. From the preenvelope of the received signal 

= ? -i0 

YAD=X (t-te (42) 
the real received signal can be described by 


y,(t)=Re{# (t-te 


(43) 
=x(t-t,)cos(@_)+ix(¢-T, )sin(@, ). 


E. THE PARABOLIC EQUATION MODEL 


The parabolic equation model used here is based on an approximation to the 
Helmholtz description of the wave equation in a cylindrical coordinate system. The majority 
of the ocean environment is well suited for a description in cylindrical coordinates. 
Assuming a time harmonic solution, the Helmholtz equation in cylindrical coordinates takes 


the form (see, for example, Jensen et al., 1994) 


18 


ic — iu SP{t2Z.) Opfrz,) 
er or rr? A? az? 


i 


+ko n *(r,.2,))p rKT,Z, b)=-4nP,O(F-F,), (44) 


where 


PUZOSD =pirz.oje™", (45) 


W . = a 
9=—— 18 a reference wave number and the acoustic index of refraction is defined by 


Co 


c(r,z,) 





(46) 


As energy is primarily radiating outward from the source, p/7,z,@) can be 


approximated by 


pfr.zo)=¥ {7,2,6)H, (kyr), (47) 


where ¥/r,z,Q) is a slowly varying envelope and H, (k,r) 1s the outward going Hankel 
function. In the far field, the Hankel function takes the form of an outward going plane wave 


and so we define 


pfrz)-— W020 “ (48) 
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Substituting this into Eq. (43) yields 





2 2 2 
ay Sat NN Sn TA a) (49 
2 8? az? Ak r2 ) 


+ 
or? "ér p2 Q kop 


where the source function on the nght-hand side has been dropped. Neglecting the azimuthal 
coupling and the near field terms and assuming that is a slowly varying function with 


range, we arrive at 


: 2 tk 
SS ae. (50) 


Defining the operators 


Hee)" 
LT =-—| -ik, — 
fog a 
and 
os = (2-1) 
2° 5 ; (52) 
this can be written as 
Pero ts 
they (Lp Uap) E. (53) 
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With the operators defined above, this constitutes what is commonly referred to as the 
“standard” parabolic equation (Tappert, 1977). For this work, we have employed the higher 


order “wide angle” parabolic equation (Thomson and Chapman, 1983) with operators 


defined by 
7 2 \ 12 ee 
(bn ene | ee) ee tee (54) 
el =| =) 
and 
Uyape= (I), (55) 


The parabolic equation model that was used for this project is the University of 
Miami Parabolic Equation model (Smith and Tappert, 1994). A split-step Fourier algorithm 
is used to numerically integrate the solution in range. This involves alternatively applying 
the U,, and the 7, operators in the z-domain and the k,-domain, respectively, where each 
operator is simply a scalar multiplier. The algorithm for stepping in range from rtor+ Ar 


can then be succintly expressed as 


P(r +Arz)=e MEET (eM CETL (7,2), (56) 
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where the wide angle 7 op Operator in the k,-space 1s defined as 


: k 2)” 
Fvapp(k,)=1- | z | (57) 


0 


The output of the model is in the form of the field functions (magnitude and phase) and 
has been referenced to a source strength of 1 uP at 1 m. The field 1s only computed at 
discrete points and the spacing in depth and range are the primary parameters that determine 
the accuracy of the results. The wavelength of interest in this study 1s 1.5-2 m and the grid 
spacing used 1s 40 cm in depth and 50 cm in range. Note that the grid size is ~ A/4 which 1s 
necessary to obtain the highest accuracy up to + 90 degrees propagation angles. In general, 
we would expect that a courser grid size would produce similar results in shallow water 


environments where propagation angles are often limited to < 15 degrees. 
F. PE ARRIVAL STRUCTURES 
The pressure 1s defined in terms of the PE field function by 


ne 0 ikor 
Pr.zf) Fe sivecaii): (58) 


where P, 1s the amplitude of the source pressure measured at R=1 m. The broadband results 


were obtained by running the model multiple times for all discrete frequencies in the 
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bandwidth under consideration and then performing a Fourier synthesis to yield travel time 
results. Assuming a normalized source amplitude, the complex arrival structure of the 


pressure field can then be written as 


P(rz.t)= | p(rz feat 


—coo 


; (59) 
= fe" War fie™*af 
r ~oo 
- i2nf— r 
Since e'” =e  “°, this phase factor can be neglected by defining the reduced time T=t-— 
“5 
such that 
l t 21 
P(rz.1)=— [ Per.oye "af (60) 
Vee 


Note that using a reduced time has no influence on the autocorrelation function. 
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il, LOCALIZATION ALGORITHMS 


The Transient Localization Project at the Naval Postgraduate School has in recent 
years studied different localization algorithms for the scenario of one receiver hydrophone 
and a point source. All the routines were based on fundamental concepts of generalized 
correlation functions described in Bendat and Piersol, 1996 and most of them are described 
in Miller et al., 1996. 

Localization algorithms may be considered generalized beamformers in which the 
plane wave replicas have been replaced by more complicated replicas of the acoustic 
propagation (e.g., modes, beams, or the vertical pressure field). The algorithms, usually 
referred to as processors, are in most cases based on a Hermitian quadratic product. The 
exact form is determined by the constraints that are put on the processor output. 

The main algorithm developed previously developed at the Naval Postgraduate 
School was the Signal Integration Filtering in Time (SIFT) algorithm which is a form of time 
domain autocorrelation matching (Benson, 1995). This algorithm had considerable success 
in the localization experiments 1n the Barents Sea (Benson, 1995). Although the results were 
favorable, there were still questions about the influence of environmental model parameters 


and acoustic (ray) model parameters on the results. 
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All the localization algorithms described below are described for the single 
hydrophone case but can easily be extended to an array of hydrophones. The source is 
described as an omnidirectional point source and the receiver is a single omnidirectional 
hydrophone. 

All localization algorithms are based on matching the received signal with a replica. 
Replicas are simulated received signals of synthetic sources at a large number of gridpoints 
in a search space. When the propagation model is 100% correct and the synthetic and real 
source functions are equal, the real and synthetic received signal should exactly match when 
the real and synthetic source positions coincide. Results of the localization are usually 
represented as an ambiguity surface. All MFP algorithms are based on this principle in some 
form. Using the reciprocity principle reduces the amount of work to a manageable size, and 
makes localization possible in a reasonable and operational feasible time. The reciprocity 
principle states that in an environment without time variations (e.g., currents) the acoustic 
pressure at location B from an omnidirectional source at location A and the acoustic 
pressure at A from an equivalent source at B are inversely proportional to the densities at 
A and B. If we assume the densities are the same at A and B, we may interchange 


hydrophone and source locations and thereby reduce the amount of work. 
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A. FULLY COHERENT LOCALIZATION 


Fully coherent localization is based on a normalized time domain cross correlation 


of a signal and the set of replicas. Given a received signal 


r(t)= ii R(fyei?™ df, (61) 
and predicted replicas 
p(t)= [Ped (62) 
the cross correlation is given by 
[R “(fP(fye laa 





fiRnPar[Ponear 


In the frequency domain, the replicas may be represented by 


PN=HPS/) (64) 
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where H(f) represents the propagation model and S(f) represents the source function of the 
replica, or template signal. 

The maximum of the peak will reflect the degree in which the modeled signal 
matches the real received signal. Only the peak value of this cross correlation is utilized. 
When the synthetic and real source coincide and the ocean model is 100% correct the value 
should be 1. This method suffers from a number of flaws. R() and P() are in general 
complex and their phase is very dependent on the ocean model parameters. (We neglect the 
source function dependence which will be discussed in the chapter on future work). A phase 
mismatch of 90 degrees at a given point will produce a zero cross-correlation. Such a 
mismatch at 1000 Hz can easily be caused by small sound speed or bathymetry errors, so this 
method does not seem to be viable for practical situations. Applying the algorithm at 


basebanded signals does not solve the phase problem. 


Be SEMI-COHERENT LOCALIZATION 


To avoid the problems with phase mismatch, which are more likely at high 
frequencies, a semi coherent approach can be taken. In this approach, the absolute value of 
the time domain replicas and the received signal is passed through a lowpass filter. This is 
similar to broadband demodulation. The matching algorithm is then performed on the low 


pass signal, 1.e., 
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R=) [| [ RO)?" “afle “2TH (65) 
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P=) {| [P07 ape at (66) 
and so 
[R “f\P “fje?™ af 
CO) (67) 
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The match is again found at the maximum of C,,(7). Previously this algonthm was used 
without much success (Miller et al. 1996) and will not be used here. However, careful 
analysis of the effects of errors in the parameters of the environment or the eigenrays on the 
atrival structures may shed some light on the poor results. 

The replica P ‘(f) shows that the largest amplitudes coincide with the earliest arrivals. 
This portion of the signal therefore provides the largest contribution to the result of the 
match. However the amplitudes of these arrivals are most prone to error due to poor 
estimation of the amplitudes of refracted rays which form the largest part of the initial 


arrivals. The semi-coherent algonthm may avoid the phase errors but amplifies the 


29 


amplitude errors. An improvement could be to normalize the amplitude of the replica. A 
number of other improvements could also be made, such as sampling at the arrival times of 


the individual rays, or clipping the template. 


C. TIME-DOMAIN AUTOCORRELATION MATCHING 


One of the problems with the previous algorithms is the lack of an absolute time 
reference. The time domain autocorrelation matching algorithm removes the absolute time 
reference from the received signal and the replica. The autocorrelation of the received signal 


is defined as 


Gaglt)= [RORHe?™ af | (68) 


and the autocorrelation of the replica is similarly 


Gpplt)= [ P“(NPNe af. (69) 
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The autocorrelation matching can then be computed by 








GGG (1) 
ise ee : (70) 
DACHO Pe Cae)? 
or in the frequency domain, using 
Ga(= [ Gag(tie "af, (71) 
Gpp(f)= f Gppt)e ?™df (72) 
then 
Arp (73) 


YD Ge DPX Geol F 


Originally the matching was designed for the real signals in the time domain. 
There are two primary advantages of autocorrelation matching. First the influence 


of noise tends to concentrate near zero lag values. By cutting out values around zero lag, the 
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effect of noise is significantly reduced. Cutting out values around zero lag also has a distinct 
effect on the contrast of the ambiguity surface as can be seen in Fig. 2. Matching values now 
range from -1 to 1. Values between -1 and 0 are neglected in all subsequent plots. The 
second advantage is the removal of any absolute time of reference. Only relative arrival 
times of separate multipaths are compared. 

Throughout most of this project, the received signals will be generated synthetically 
using the PE model. Replicas were generated using the same PE model for various 


environments to determine the optimum number of lag values near zero to remove. The 





Figure 2 Effect of notching out near zero lag values of autocorrelation matching at a 
certain range. Plots (a) - (f) show increasing number of notched out values. 


a2 


results indicate that this number depends on the shape of the source signal and therefore on 
the bandwidth. Furthermore, where it has a positive effect on the contrast, it has a negative 
effect on the footprint size. For example, by removing seven points near zero lag (zero lag 
plus three at positive and three at negative lags) for a negative gradient environment, the 
footprint is reduced from about five meters in depth to about two and a half meters in depth. 
The footprint size does not seem to be reduced further when more points than the optimum 
number of points are removed. However, the removal of more points begins to degrade the 
performance by increasing the level of false peaks (sidelobes). Note that the removal of near 
zero lag points in the time domain is equivalent to removing the mean and lowest order 
trends from the frequency domain. We found that the removal of points corresponding to a 
width of about 2/bandwidth is adequate under most circumstances for source signals having 


a Blackman shaped frequency response. 
D. SEMI-COHERENT AUTOCORRELATION MATCHING 
As has been done for semi-coherent localization a less coherent form of localization 


can be defined for the autocorrelation matching. Using Eq. 71 and Eq. 72 we can define a 


semi-coherent autocorrelation matching algorithm as 
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This form reduces the influence of the phase. Working with only positive values, however, 
makes the dynamic range of the ambiguity surface very small. Results for PE generated 
synthetic received signals and PE generated replicas look promising but the performance 
when comparing PE received signals and ray model replicas was not adequate to pursue this 


algorithm at this time. 


E. THE SIFT LOCALIZATION ALGORITHM 


The SIFT (Signal Integration Filtering in Time) algorithm has been developed as a 
result of previous research (Miller et al. 1996, Benson 1995). It is based on a ray code 
solution. The templates are generated using the arrival time and phase of each eigenray. 
These parameters determine uniquely each gridpoint in a search grid. The SIFT algorithm 
assumes that the signal can be written as the sum of scaled, phase and time shifted versions 


of the original signal, 1.e., 
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N 
WO=d, y,(t) 


5 Re(a x tot Je. . °n) (75) 


55 a,[x(¢-t, )cos(8,)+x(¢-T, )sin(O_)], 


where a, 1s the real amplitude and t, and 0, are the travel time and phase, respectively, of 


elgenray, x(t) is the real source signal and x(¢)is its Hilbert transform. The autocorrelation 


of the replicas is given by 


G.y(2)= [ Oye) 


or (76) 
DEP 
n=1 m=) a 
The correlation between individual eigenrays can be wnitten as 
Gy y (T= [a Vm(E*T) (77) 


Or 


(t)=a a _ cos® cos® R_(t+t, -T,) 

: +a.a_cos@ sind R(t+T, -T,,) 
+aa sin8 cosO_ R.(t+T,-T,) (78) 
+a.a_sin® sin8 R.(t+T, -T,). 
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Using properties of the Hilbert transform, this can be rewritten 


taa G_(t+t,-t_)sin(8 -6), i) 


nm om XX 


G, . (t)=a,a_ G_(t+t,-t, )cos(8, -9,) 


where G.. is the Hilbert transform of the real source autocorrelation, G,, 
In the original SIFT algorithm, three major approximations are made: 
1) the amplitudes are set equal to unity (a,=1 for all n); 
2) the phases are independent of frequency; 


3) the autocorrelation functions, G. and G are approximated by discrete 
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functions defined by 


CAG face S84) + . 15 le ~T (80) 


otherwise 


and 


eta te. aeOt 
CA GT Te) el T=T_—T 407 (81) 
0 otherwise 


where Ot is chosen to be the correlation time of the measured signal. 
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The matching part of the algorithm 1s then defined by 


GaG 
= PP RR 
(82) 


/IG,,I7IG, 1° 


where R,, and Rep are the vectors of the sampled correlations of the replicas and the received 
signal 
-a,a,sin(9, -9,) 
a,a,cos(O, -9,) 
a,a,sin(0, -9,) 
Cai -a,a,sin(0, -9,) (83) 
a,a,cos(8, -9,) 


aa,sin(O, -9,) 


a 


L 


and 


Gate ts Oy 
Gpg(T -T2) 
Get, -t,+5T) 
Ce GEG -7,-Om)e (84) 
Gpp(t, -75) 
Gaal, —t, +00) 


L : | 
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Usually the autocorrelation of the received signals will be formed as a linear correlation and 
will then be sampled. The above description assumes the eigenrays are sorted in order of 
increasing travel time. 

A similar derivation can be made for the representation of a complex signal. The ray 


arrival structure can then be represented by 


N 
VWO=Y FO 
n=] 
N = f 
= a [R(t-t Je Mn ny (85) 
n=] 


. | 
=)> a [X(t-t, cos(2nf.t, +0 )-ix(t-1,)sin(2nf,t, +8 ) 
n=] 


In general, we may want to use signals that have been basebanded or shifted in frequency 


by some amount /,. The complex source may then be defined by replacing 


R(t-t1)X(t-t)e 77 (86) 
The autocorrelation of the received signal can now be expressed as 


G5= (HOF) 
NON 


=e G, (0) 


n=l m=1 


(87) 


38 


To simplify the expression we define 


@,=207,, +8, (88) 
and then the correlation between the signals of the individual eigenrays can be expanded as 


Gi 5 (t)= [¥,Oyg(t+a)alt 
= aa, cos(h)cos(, Ret tia ca) (89) 
- aa icos(P,)sin(®d, )R.{t+T -T) 


+ aa_isin(®, )cos(p, )R{t+T,-T) 
+ aa sin( )sin(® )R.{t+T -T) 


Or 


G; (t)=a,a G.{t+t -t, )[cos(®, - )-isin(d, -,)] 
= apa (Gaelo te. Co LCOS 2 70ha( Cire) dae a) (90) 
-isin(2ny (t-t )+0_ -@ Df. 
If we would apply similar approximations as with the original SIFT algorithm we would 
approximate the amplitude with unity and the autocorrelation function would have to be 
sampled with a single point approximation for t=T,,-t,. When taking only the positive or 


only the negative lags into account for the approximation, the matching function would be 


GAG 


S=Re ~ PP RR 


{1G PIG. I 


(91) 
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When both the negative and positive lags are used the result will obviously be real. For 
bandpass signals this seems a worthwhile approach. An alternative approach to 
accommodate bandpass signals would be to shift them to the lowest possible frequency 
band, form the real autocorrelation of the received signal and then apply the onginal SIFT 


algorithm for real signals. 


F. SIFT AND THE PHYSICS OF THE LOCALIZATION PROBLEM 


The choice of dt is one of the major difficulties with the SIFT algonthm. One way 
to optimize the algorithm would be to try and optimize the results for a small range of 6t’s. 
For simulations dt can easily be calculated from the cross correlation of x(t) and x(t). The 
performance may be enhanced by optimizing over a small range of 6t’s, although this 
reduces the computational speed. For this approach to be valid, it is also necessary that of 
is less than the separation between distinct arrivals and greater than the time sampling rate. 
For very broadband signals with good time resolution, this approximation may be very good. 
However, this method is quite crude and may not work for narrow band signals. An 
alternative may be to use a multi-point approximation instead of a three-point 
approximation. 

The 3-point approximation of the correlation is a low frequency approximation. At 
high frequencies and for bandpass signals this approximation is no longer valid. Furthermore 


the values with which the approximation is made (1,-1,1) are dependent on the shape of the 
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source function. It may be worthwhile to look at a less ngid approximation. A few test cases 
were examined in which the Hilbert canons of the autocorrelation was neglected and no 
serious degradation in performance was observed. 

The phase is one of the most sensitive parameters used in the algorithm. Small 
changes in phase at turning points and boundary interactions give the phase at larger 
distances a more or less random value. Although arrival structures may look similar, the 
autocorrelation will be very different, and much of the match is lost. For the complex 
envelope version of SIFT the result not only depends on the phase difference due to 
reflections and refractions but also on the phase difference due to the time delay. Although 
the travel times are known to be a very robust physical parameter in these cases, small 
errors may become significant at high frequencies. 

In general, the calculation of eigenray amplitudes is reasonably accurate. However 
there are a few arrivals in those environments where the estimation of amplitudes is quite 
poor. Typically, these rays correspond to refracting rays near caustics producing erroneously 
large amplitudes. A large arrival will dominate the autocorrelation and make it look very 
much like the arrival structure. This may be a good reason to normalize the amplitude. When 
you normalize the amplitude, the later arrivals have the same weight in the correlation and, 
therefore, in the matching algorithm as the earlier arrivals. This has a number of advantages. 
The later arrivals carry as much information as the first few and, therefore, taking them into 
account increases the amount of information in your match. Also solutions of the eigenray 


extraction routine are more stable for the later arrivals than for earlier arrivals. 
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IV. EXPERIMENTAL SETUP AND RESULTS 


To gain insight into the performance of the autocorrelation matching algonthms and 
approximate autocorrelation matching algorithms, several computer experiments have been 
done. The experiments have primarily been performed for a frequency band of 750-1000 Hz. 
The baseline for the results was defined as the results of the autocorrelation matching of a 
PE generated source and PE generated templates. For most of the experiments a source 
range of 5500 m, a source depth of 59.8 m and a receiver depth of 40.2 m have been 
assumed. For the environments defined below, at this range the first arrivals arrive at 
approximately 3.7 seconds after the transmission and the significant part of the arrival 
structure is about 0.4 seconds in length. As the project did not focus on the frequency 
characteristics of the source function, the same shape has been assumed for the source signal 
and the templates. Throughout the project analysis the absolute scale of the pressure field 
calculation has been neglected as the similarity of the arrival structures 1s emphasized not 


the absolute values of the pressure field. 


A. ENVIRONMENTS 


For the synthetic experiments, three shallow water, range independent environments 


have been defined. The first has a negative sound speed gradient (1500-1475 m/s), the 


second has aconstant sound speed of 1500 m/s, while the third has a positive sound speed 
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gradient (1500-1501.5 m/s). All other acoustic parameters are fixed for the three 
environments. The water depth is 100 m and the bottom parameters, the attenuation, 
density and sound speed, were chosen as «=0.48 dB/(km Hz), Pronom=1900 kg/m’, 
Croton = 1650 m/s, respectively, consistent with a sand bottom. The bottom 1s modeled as an 
infinite half space in both propagation models. The propagation loss for the environments 
at a frequency of 875 Hz is shown in Fig. 3 (a), (b), and (c). These plots suggest a very rapid 
variation of the acoustic pressure with respect to depth. The vertical size of localization 
footprints is therefore anticipated to be small. 

The third, positive gradient, environment has only been used to generate a limited 
set of ACM baseline results. It was, perhaps wrongly, anticipated that the results might not 
significantly differ from the results generated for the iso-speed environment, and therefore 
the computational effort would not outweigh the additional results. While this appeared to 
be the case for the third environment (with a weak positive gradient), the results for the first 
environment (strong negative gradient) were significantly different. The analysis suggested 
that most of the algorithms performed well for the second environment (iso-speed) but were 
at their performance limits for the first environment. It is therefore likely that a less severe 
sound speed gradient environment could have given more insight into the exact performance 


limits of the algorithms. 
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Figure 3 Transmission loss at 875 Hz for (a) a negative 
sound speed gradient (1500-1475 m/s) environment, (b) an iso- 
speed (1500 m/s) environment, and (c) a positive sound speed 
gradient (1500-1501.5 m/s) environment. 
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B. HAMILTONIAN RAYTRACING PROGRAM FOR THE OCEAN 


The Hamiltonian Raytracing Program for the Ocean (HARPO) has a number of 
parameters which directly influence the accuracy of the results. To be able to use the results 
for this application, a maximum ray launch angle separation of 0.05 degrees has been 
determined. Larger spread caused unacceptable inaccuracies in the eigenray extraction 
program. The program has been set to output all the points, and the integration accuracy has 
been set to 10°. Less accurate settings would degrade the performance of the eigenray 
extraction program. Travel times and eigenray calculations have been compared with results 
from an analytical model and proved to be very accurate. To reduce the computational 
burden, the assumption has been made that at the ranges of interest the attenuation due to 
bottom reflections would be high enough such that source angles larger than critical could 
be neglected. This reduces the number of eigenrays to about 40-50 per grid point for the 
negative gradient environment for a range of about 5 km. For shorter ranges, source angle 
limits of +40 degrees have been used resulting in more than 18 eigenrays per grid point for 
ranges above 1 km. This set up resulted in a ratio of the smallest to the largest eigenray 
amplitude per grid point of smaller than 0.02 for 95% of the environment. The grid point 
spacing of the eigenray extraction routine has been set to 5 m in range and 1.4 m in depth. 


The depth spacing does not meet the initial minimum requirement of five grid points per 
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dimension of the footprint (12.5 m x 1.2 m at 5.5 km for a source depth of 59.8 m and a 
receiver depth of 40.2 m). It is a compromise due to the computational speed of the 
eigenrays extraction routine. 

During the experiments, four eigenray extraction problem areas were indicated. Very 
abrupt change, or little change, together with noise due to round off and interpolation errors 
and using too large angle separations creates either too few or too many eigenrays or 
‘spurious’ eigenrays at these points. For the iso-velocity environment, a separate program has 
been used that is based on an analytical solution. This allowed a little more flexibility and 


considerably reduced the run time. 
C. UNIVERSITY OF MIAMI PARABOLIC EQUATION MODEL 


The University of Miami Parabolic Equation Model (UMPE) was used with the 
wide-angle PE approximation, as described in Chap. I. The maximum computational depth 
has been set to 409.60 m. A depth spacing of 40 cm leads to 250 grid points within the 
water column, and 4-5 grid points per wavelength. As a comparison, only 70 grid points in 
depth have been used to cover the water column for the ray model. The grid spacing in range 
has been set to 50 cm which is about the minimum for an accurate prediction in the 
environments described above. Results have been stored in 2.5 m and 5 m range increments 
for the negative gradient and iso-speed/positive gradient environments, repectively. The 


frequency is sampled from 750 Hz to 1000 Hz in 257 frequency bins. At a range of 6500 m, 
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this frequency resolution still allows the calculation of the correlation as the Fourier 


transform of the power spectral density without degradation due to wrap-around. 


D. RECIPROCITY 


In the experiments, the source and receiver are assumed to be a point source and a 
point receiver. The source spectrum is shaped as a Blackman window. This avoids large 
sidelobes in the time domain arrival structure and leads to a pulse length of approximately 
8 ms. To optimize the eigenray extraction routines, a large number of reciprocity checks for 
different source and receiver locations have been done. The autocorrelation function of the 
arrival structures proved to be an excellent tool in cases where the arrival structures looked 
alike. Small differences in amplitude and especially phase led to large differences in the 
autocorrelation function. This suggested that localization algorithms based on the time 
domain autocorrelation would be very sensitive to these parameters, particularly the phase. 

The PE-results showed very good reciprocity. The ray code results showed good 
reciprocity except at certain locations where erroneous amplitude calculations degraded part 
of the arrival structure. Looking at the difference between the ray code and PE arrival 
structures, it appears that the largest difference can be found in the first few arrivals. The 
tails of the arrivals structures are very similar. This suggests that the ray code approximation 
for refracting rays is less accurate than for reflecting rays. Perturbation analysis of the phase 


and amplitude of the rays suggests that the arrival structures and autocorrelation are more 
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susceptible to phase errors than to amplitude errors. Results also suggest that certain parts 
of the autocorrelation are more affected than others. Very short lags and very long lags seem 
to suffer more than intermediate values. Knowing that short lag values are notched out and 
large lag values will have no significant contribution to the match suggests that the 
algorithms should at least have some degree of tolerance toward amplitude and phase errors. 
This suggestion is based on visual analysis of a large number of results, but no statistical 


analysis has been performed. 


E. MATCHING ALGORITHMS 


Numerical experiments have focused on autocorrelation matching and the SIFT 
matching algorithm. The baseline for the results has been defined by the PE autocorrelation 
matching over the full 250 Hz bandwidth using a Blackman shaped frequency response. The 
autocorrelation matching experiments can be separated into their variations of the basic 
algorithm. The following variations have been specifically examined: 1) the effect of 
notching out points near zero lag; 2) the effect of using a smaller bandwidth, which was 
anticipated to give a larger footprint size; 3) the effect of splitting up the frequency band into 
several smaller frequency bands; 4) the effect of a lower frequency resolution; 5) the effect 
of raising the tail of the autocorrelation function, i.e. to increase the weight of large lag 


values (corresponding to cross terms of early and late arrivals); 6) the effect of matching the 
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absolute values of the autocorrelations functions; and 7) the effect of amplitude and phase 
perturbations on the frequency domain arrival structure. 

Specifically for the ray results, additional effects have been examined: 1) the effect 
of raising the tail of the arrival structure to increase to weight of the later arrivals; 2) the 
effect of neglecting the first few arrivals, as the amplitude and phase of these arrivals are 
suspected to be the cause of most of the bad results; and 3) the effect of matching the arrival 
structures directly, in both complex and absolute form. 

The SIFT algorithm has been considered in four different ways: 1) a form where the 
complex signals have been base-banded and a 1-point — to the autocorrelation 
function is used; 2) a form where the real signals are shifted in frequency to a center 
frequency of half the bandwidth and where subsequently a modified form of the original 
SIFT algorithm is employed; 3) a form where the signals are shifted and then a zero mean 
version of the original SIFT algorithm is employed which results in a multiple-point match; 
4) variations on the original low frequency implementation of the SIFT algorithm. For all 
the variations, notching out points near zero lag, neglecting early arrivals, and variations of 
the parameters that describe the approximation to the correlation functions applied in SIFT 
have been considered. 

The grid spacing of the templates has been based on the footpmnt size of the baseline 
results and computational arguments. As the main area of interest is beyond 5 km, the results 
of autocorrelation matching at 5.5 km for the negative gradient environment have been used 


for determining the footprint. The footprint size 1s approximately 1.2 m in depth and 12.5 


50 


m inrange. With a norm of 5 points per dimension of the footprint this would lead to a gnd 
spacing of 20 cm in depth and 2 m in range. Analysis at other ranges suggests the footprint 
will be smaller at smaller ranges and larger at larger ranges. The excessive run time of the 
eigenray extraction, the available disc storage and run time for the matching aigorithms have 
iead to a grid spacing of 2.5 m in range and 40 cm in depth for PE generated templates in 
the negative gradient environment, 5 m in range and 0.4 m in depth for PE generated 
templates in the iso-speed and positive gradient environments, and 5 m in range and 1.4 m 
in depth for the ray code generated templates. For the negative gradient environment, the 
eigenray results have only been calculated for 1000-2000 m and 4500-6500 m due to the 
computational load of the MATLAB based eigenray extraction routines. Attempts to apply 
the recently available MATLAB c-compiler were not successful. The presently available 
compiler still uses many standard MATLAB routines which prevent generation of efficient 
code. Also the nature of the problem which relies on using complex numbers degrades the 
ability to generate efficient code. 

To generate synthetic ‘received’ signals, three basic algorithms have been used: an 
impulse model based on Eq. 43, a time domain model based on Eq. 33, and a frequency 
domain model based on Eq. 37. Where the impulse model is an inherently low frequency 
modei, the time domain model is suited primarily for band pass signals. A small 


modification to the original time domain based model developed by Chiu made it suitable 


for low frequency applications. In the ttme domain model, the tails of the pulse are removed 
which is a major advantage in testing algorithms where overlap of pulses is a crucial 
parameter. 

The autocorrelation of the synthetically generated received signals has been 
computed using several different methods. For autocorrelation matching using frequency 
based data as generated by the UMPE model, it may seem that the transform of the power 
spectral density may be the most practical solution. However the circular convolution does 
not work well with pulse signals where the pulse length is more than half the transform time 
base. For this reason the autocorrelation for this project has primarily been computed from 
the time domain arrival structure. Where in previous work an ‘unbiased’ autocorrelation 
estimate has been used, for this work the ‘biased' autocorrelation estimate has been used 
which is more suited for pulse shaped signals. As the time base is determined by the 
frequency resolution, it may be beneficial in future experiments to artificially stretch the 
time base before performing the autocorrelation. Results of applying phase disturbances to 
the PE results suggested that keeping the ttme domain arrival structure centered in the first 
half of the time base can not always be accomplished. While theoretically PE phase errors 
should have no influence on the autocorrelation results, there will be a degradation in the 
results due to the above. Thus, while theoretically the time domain based autocorrelation 
may be preferable, in practice the frequency domain based calculations are more robust. 

No use has been made of the usual performance measures such as peak to sidelobe 


level or mismatch in range and depth. The reason is that these measures tend to put a 
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number to a match. This work tries to explain how and why the match will degrade. It has 
been observed that a match tends to become amplified along certain lines of high 
correlation/energy. Furthermore, when the matching parameters exceed certain values, 
depending on the type of matching, environment, range and bandwidth, the results suggest 
the single match will explode into several 'matches' spread around the actual source location. 


This still gives an indication of the source location but it will not give a single point match. 


F. COHERENT MATCHING RESULTS 


Coherent matching has been performed in the time domain and in the frequency 
domain. The frequency domain approach may require an additional parameter to allow for 
time shifts. All presented analysis uses a time domain approach. Results may slightly differ 
due to the time resolution used. Results for the negative gradient environment with a source 
at 5500 m display sidelobe levels comparable to autocorrelation matching with an optimum 
number of lag values notched out. The footprint size for the same setup was found to be 
larger for the coherent matching than for the autocorrelation matching. Coherent matching 
results using ray code templates and PE generated source signals showed some phase shift 
between the templates and the source signals. In general, we were able to localize correctly 
up to ranges of more than 5 km. Short ranges resulted in very high sidelobes, especially in 
the 1so-speed environment. No significant changes have been made to the basic coherent 


matching algorithm. In many cases where the autocorrelation matching algorithms failed to 
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localize, the coherent matching still gave reasonable results. Taking into account a phase 
shift in the result of the match improved the ne significantly. There are however 
many drawbacks to coherent matching, and there are locations where we found that many 
templates and the received signal were almost alike which led to very high sidelobes. 
Discretized phases may be a convenient way to improve the coherent matching algonthm. 

In addition some analysis has been done on a less coherent approach using the 
absolute value of the time domain arrival structure. The dynamic range was improved by 
subtracting the average over a certain area of the matching resuls. Altough baseline results 
showed promising, the results using PE generated received signals and ray code generated 
templates showed very high sidelobes which prohibited correct localization. A similar 


approach will be discussed further in Chapter V. 
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Figure 4 Coherent matching for the 
negative gradient environment 
(baseline result), with a source at 
S500 m range and 59-6 m Geeth ancme 
receiver at 40.2 m depth. 
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Figure 5 Coherent matching for the 
negative gradient environment using 
a PE generated source signal and ray 
code templates with a source at 5500 
m range and 59.8 m depth anda 
receiver. at 40.2 Mm Gaocn, 2eceuneiig 
for the phase shift results in 
correct localization. 
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G. AUTOCORRELATION MATCHING RESULTS (BASELINE) 


The autocorrelation matching algorithm is described in Chap. III, Section C. To 
improve upon the results of the basic algorithm, the following variations have been 
considered: notching out near zero lag values; enhancing the dynamic range by subtracting 
the average over a certain area of the match and renormalizing; increasing the significance 
of large lage value by filtering; changing the sample. rate, bandwidth and frequency 
resolution; and splitting the frequency band into several small bands. To get an initial 
assessment of the robustness, amplitude and phase perturbations have been added to the PE 
generated received signals and PE generated templates. 

Results for both the negative and positive gradient environments and a source at 
5500 m range suggest that notching out zero lag and the first 2-4 lag values (corresponding 
to (1 or 2)/BW) gives an optimal performance in terms of peak to sidelobe level. For the 
| isospeed environment the optimum was found by notching out four times the previous 
number of lag values, having still larger sidelobes than the other two environments. 
Presumably arrival structures are more similar for the iso-speed environment. Results 
suggest that the optimum number of notched out lag values depends on range and, to a lesser 
extent, on bandwidth (except for low frequencies where the bandwidth dominates). The first 
few lag samples represent the low frequency trends in the autocorrelation. They are more 
related to the power in the signal rather than the shape and will dominate the match when 


they are not removed. Notching out these lag values, however, decreases the footprint size. 
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Notching out more than the optimum number of lags does not significantly reduce the 
footprint size any further for the negative gradient environment. For the iso-speed 
environment, however, the footprint size continues to decrease further when the number of 
notched out lag values increases. 

As an alternative to notching out values around zero lag in the time domain, two 
approaches have been considered. Removing the mean from the frequency domain results 
leads to similar results as notching out only the zero lag value. Removing lower order trends 
would require filtering of the power spectrum, which makes it much more complicated than 
notching out values in the time domain, with theoretically the same results. The second 
approach aims to increase the dynamic range, and is based on the subtraction of the 
matching value, averaged over an area, from the individual matching values. Contrary to 
notching out values around zero lag, this approach has no real physical interpretation. The 
results are comparable. However sidelobe levels are slightly higher, and it does not have the 
degrees of freedom to optimize the results as with notching out lag values in the time 
domain. 

The footprint was found to depend on range, receiver/source depth, sound speed 
profile, position of the source relative to the ‘lines of high correlation/energy', and on 
bandwidth. It was found for the cases considered that the general trend is for the footprint 
to increase with range. Footprints for the same source/receiver depth were found to be much 
larger than for significantly separated source receiver depths. In general, footprints tended 


to be larger for the environment with a small positive gradient and the iso-speed 
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environment than for the environment with a large negative gradient. The position relative 
to a line of high correlation/energy determines if and how the footprint will stretch. The 
most significant difference between the slight positive gradient or iso-speed environment on 
one hand and the severe negative gradient environment on the other is the vertical extent of 
the footprint which is nearly twice as large for the slight positive gradient and even larger 
for the iso-speed environment than for the negative gradient case. The smaller footprints of 
the negative gradient eeocrnent can be directly related to the more severe interference 
patterns. The significantly larger footprints for the same source receiver depth are thought 
to be caused more by energy matching than by matching of the arrival structure shape. 
The number of points or sample rate did not appear to be a crucial parameter. There 
are some results suggest that a sample rate of 1000 Hz might give better results and a slightly 
larger footprint than a sample rate of 250 Hz. For this reason the bulk of the analysis has 
been performed for a sample rate of 250 Hz for the negative gradient environment and 1000 
Hz sample rate for the iso-speed environment. Contrary to our initial expectations, smaller 
bandwidths generally produced smaller footprints and higher sidelobes than larger 
bandwidths. This may be due to fewer arrivals being sampled at smaller bandwidths cause 
the the autocorrelation (and its match) to drop off rapidly away from the source location. 
Larger bandwidths on the other hand, resolve more multipath arrivals which may broaden 


the autocorrelation and generate a larger footprint. 
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A lower frequency resolution showed a larger footprint but also large sidelobe 
levels. Although the number of cases examined 1s limited, the results suggest that larger 
bandwidths are preferable, as expected. Results also suggest that it may not be necessary to 
apply any window. Results show that the degradation due to sidelobes/ringing in the time 
domain may be outweighed by the increase in performance due to the larger bandwidth. 

Separating the 750-1000 Hz band into 4 smaller bands and averaging the results led 
to a slightly larger footprint as compared to the full bandwidth results but with a significant 
increase in sidelobe levels. Later arrivals are more separated in time than earlier arrivals and 
are also more stable in these environments as they are not related to refracted propagation 
paths. To be able to exploit them, the later arrivals should be enhanced. In the PE model the 
individual arrivals cannot be manipulated. Therefore, an alternative approach was used, 
enhancing the tail of the autocorrelation. Results suggest that indeed some improvement may 
be gained with this approach in terms of the sidelobe levels, especially if it is applied to both 
the received signals and the templates. The footprint, however, appears to be smaller than 
for the full bandwidth case (7.5 m in range and 2.5 m in diameter at 5.5 km for the negative 
gradient environment). The results were not found to be significant enough to pursue this 
approach any further. It may be that in more complicated environments this approach should 


be reconsidered. 
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Figure 6 Baseline autocorrelation matching results for the 
negative gradient environment. The source is located at 5500 m 
range and 59.8 m depth and the receiver at 40.2 m depth. Three 
notched out near zero lag values have been used. (a) Full 
range. (b) Range expanded near the source location. 
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rare 7 Baseline autocorrelation matching results for the 
1so-speed environment. The source is located at 5500 m range 
ema 59.8 m depth and the receiver at 40.26m depth. Three 
notched out near zero lag values have been used. (a) Full 
range. (b) Range expanded near the source location. 
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Figure 8 Baseline autocorrelation matching results for the 
positive sound speed gradient environment. The source is 
located at 5000 m range and 40.2 m depth and the receiver at 
59.8 m depth. No window has been applied and 3 near zero lag 
values have been notched. (a) Full range. (6b) Range expanded 
near the source location. 





To get an assessment of the robustness of the algorithms amplitude and phase 
perturbations have been applied. The effect of phase perturbations has already been 
discussed. The amplitude perturbation has been applied as a uniform random fluctuation. 
Using a random distributed amplitude fluctuation up to +/- 40 % only showed a reduction 
in the peak to sidelobe level, but no significant increase in sidelobe levels or change in the 
footprint size were observed. This indicated that the algorithm is fairly robust for amplitude 


fluctuations. 


H. AUTOCORRELATION MATCHING (NUMERICAL RESULTS) 


When matching PE source signals and ray code templates for the negative gradient 
environment using a windowed frequency response, results appear to be changing 
periodically for ranges beyond about 5 km. Moving the source over a range increment of 
4500-6500 m, a correct localization was obtained about 50% of the time. Of the remaining 
half, some regions, did not show any localization while in other places the localization was 
obscured by high sidelobes. Short ranges give more consistent results but also contain 
inherently high sidelobes. Even with a severely inadequate ray set, localization could be 


performed at ranges up to 2000 m. 
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Several changes have been made to obtain more consistent results at longer ranges. Notching 
out a large number of lag values oe improves the results but not consistently. 
Applying a window that amplifies the larger lag values of the correlation gave some 
improvement when applied to both the received signal and the templates. 

From the good localizations, a 1.2% difference in range between the ray code 
solution and the PE solution was found for the negative gradient environment. This 
difference is probably due to phase errors in the PE approximation or a sub-optimal choice 
for the reference sound speed and using below optimal values for the range and depth grid 
spacing. The footprint size 1s usually very small, 5-15 m in range and 1.4 m (1 grid point) 
in depth. Neglecting the initial arrivals in the ray code solution did not improve the results. 
It also seems that higher sampling rates give slightly higher matching peak values but do not 
significantly improve the results. Splitting the frequency band into several smaller bands can 
improve the results for a certain set of matching parameters, but the results are not 
predictable. No solid reason could be found explaining the changing behavior with range. 
It 1s speculated that this is not directly related to the bandwidth (smaller bandwidths 
sometimes produced better results than larger bandwidths) but may be related to the modal 
interference pattern, the phase errors in the initial arrivals, and possibly under-sampling. 

No correlation was found between eigenray amplitude problem areas, eigenray 
extraction problem areas, arrival times, arrival time differences, depth of the perigee of 
refracted rays, or number of refracted rays and the changing behavior with range. The phase 


of the ‘half’ match (using only positive or negative lag values) was observed to vary 


64 


periodically. But as it is composed of many parts, it is more indicative than a real measure, 
and it shows no sudden changes for locations where a match was expected but not found. 
A fixed phase difference in the time domain should not show in the autocorrelation. Still, 
as other methods do not seem to suffer as much, the reason is presumably related to the 
autocorrelation. 

The only link to an explanation may be found from results which were obtained 
without a frequency window. Using the full 250 Hz bandwidth without frequency windowing 
led to correct localizations where no localization was found with windowing. This suggests 
that any degradation due to sidelobes of the window Pessiad is outweighed by the 
improvement due to the effective larger bandwidth. An additional improvement was 
achieved by accounting for the phase shift in the ‘half match’, which led to a considerable 
enhancement of the peak to sidelobe ratio. These results may also suggest that the required 
bandwidth is not a monotone increasing function of range and the sound speed profile 
probably has a distinct influence on that relationship. 

For the iso-speed environment, good localization is found over the whole range from 
900-7500 m. At larger ranges the performance is better than at short ranges due to the higher 
sidelobes at short ranges. Footprints are slightly smaller than for the baseline results. These 
results were not found to be very sensitive to the time sampling rate, in contrast to the 


baseline results. 
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Figure 9 Autocorrelation matching using a PE generated source 
Signal and ray code templates.(a) Negative gradient 
environment with a source located at 5500 m range and 59.8 m 
depth and a receiver at 40.2 m depth. Full 250 Hz bandwidth 
has been used without window and the phase difference of the 
autocorrelations has been accounted for.(b) Iso-speed 
environment with source at 5000 m range and 59.8 m depth and 
receiver at 40.2 m depth. 
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L. SIFT MATCHING RESULTS (LOW FREQUENCY) 


The SIFT matching algorithm that was developed by Miller and others (Miller et al, 
1995), has been modified to account for bandpass signals, to bring it closer to a true 
autocorrelation matching, to improve the notching for short time lags, to adapt it for 1- to 
n-point matching instead of the fixed 3-point matching, to manipulate the amplitude of the 
cross correlation of the source function and the Hilbert transform, and to include the ray 
amplitude in the matching. The original SIFT algorithm is based on a source signal that 
behaves like a low frequency pulse signal. This type of signal has been proven adequate for 
modeling of SUS charges. However, it may be less suitable to describe other signals, 
especially signals which have been high-pass filtered. For this kind of signal an n-point 
match, where n increases with decreasing bandwidth, may be more appropriate although the 
matching may become too complicated for large n. In certain favorable cases, tuning of the 
Ot factor in Eq. 81 may even allow the use of a 3-point match when using a signal that does 
not resemble the assumed signal model. 

The original SIFT algorithm will never give a perfect match even if the received 
signals can be modeled as spikes. This is because the cross correlations are stacked as an 
array rather than summing the cross correlation of individual arrivals, . However, this result 
is matched with the sampled autocorrelation of the received signal in which the individual 
atrivals are inherently added coherently. As this addition slows down the algorithm, it has 


not always been used. 
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The notching previously was based on the arrival time difference. It did not take into 
account the approximated cross correlation of two arrivals containing two points that were 
shifted by St. In the algorithms used here this slight modification, which may be important 
at short lags, was taken into account. 

Algorithms for the 3-point matching have been implemented. The n-point matching 
algorithm has been considered but given the time scale it was not implemented. As the 
source function of the received signal was known, the amplitude of the 3-point 
approximation of the cross correlation has been tuned to optimize the performance. For an 
n-point SIFT algorithm, this would be an immense task and not conducive to working with 
real data. 

The amplitude of the rays has been brought back into the SIFT algorithms in contrast 
to previous versions. In general, this improved the results but only slightly. A few cases even 
showed a degradation of the results when ray amplitudes were included in the algonthm. 

Aside from the parameters mentioned above, the following parameters were taken 
into consideration to optimize the results: the offset 6t, the sampling rate, and the bandwidth. 

To develop some intuition about the general performance, the low frequency version 
of the algorithm has been applied to the negative sound speed and iso-speed environments. 
The results imply that the performance of the algorithm for a certain environment, source 
function, bandwidth, and bottom depth is range limited. An empirical analysis suggests that 


the area where the algorithm performs succesfully is below the curve 
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oa (92) 


where d is the water depth, c is an average sound speed, B is the true bandwidth, r is the 
range and k is a constant that depends mainly on the environment. The constant k ranges 
from approximately 0.35 to approximately 0.5 for the iso-speed and negative gradient 
environments. The approximate curves for the iso-speed and negative gradient environments 
are depicted in Fig. 10. There is strong evidence that the curves describe the areas where 
the algorithm is likely to perform, but not enough analysis has been done to exactly localize 
the curves. The formula is based on a first order approximation of the arrival time 
differences for r>>d. 

Results suggest that there may be different regimes which govern the matching. First, 
there are the high bandwidth cases which are well below the curve described by the 
parameters. It is suggested that the match consists of different contributions, most probably 
energy and shape of the arrival structure, which are influenced by the number of notched out 
values. For high bandwidths, the shape of the arrival structure prevails. When the bandwidth 
is reduced, the shape of the match stretches along the lines of high correlation/energy (most 
probably corresponding to rays with dominating amplitudes). Reducing the bandwidth 
further causes the match to split up into several ‘matching points’ which migrate along the 


lines of high correlation/energy. The improvements mentioned earlier give some relief but 
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only marginally increase the performance. Plots of the number of sampled cross correlations 
with lag values within one pulse width do not show any abrupt changes with decrease in 


bandwidth. 





range (km) (b) 


Figure 10 Estimated performance limitations for the LF SIFT algoritm for (a) the 
negative gradient environment and (b) the iso-speed environment. 


The autocorrelation of the received signal however gives a clue to the process... 
During the reduction of the bandwidth, a low frequency trend starts to dominate the arrival 


structure of the received signal. The individual arrivals seem to ride on top of this low 
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frequency trend. In the autocorrelation, the sequence of the positive/negative going parts of 
the arrival structure becomes dominant. This makes it necessary to notch out about twice the 
number of points notched out at higher bandwidths. At these low bandwidths, performance 
is not impressive and the footprints are large and only approximate localization is possible, 
but the sidelobes are excessively large. It is still unclear whether we could exploit this last 
regime, especially in more complicated or more noisy environments. The change in the 
autocorrelation of the received signal indicates the limit up to where the algorithm can be 
pushed in the ordinary form. Results for both environments considered differ but the trend 
agrees with the theory. As the significant arrivals for the negative gradient environment are 
more confined in time than for the iso-speed environment, the limits for a certain algorithm 
in this environment are more severe. No low frequency synthetic PE source signals have 
been used because results of the ray model and UMPE model are largely different below 10 
Hz due to the breakdown of ray theory at these frequencies. 

The size of the footprint is found to be dependant on range and bandwidth. In 
general, a larger bandwidth will give a smaller footprint. The numerical analysis presented 
here indicated that a bandwidth of 250 Hz at a range of 5 and 7.5 km was not adequate for 
the algorithm to succesfully localize for the negative sound speed gradient environment. 
Consequently, some of the goals of the project based on these parameters had to be dropped. 

Results show only moderate sensitivity to the offset parameter ot, and in some cases 
a forced 1-point match (neglecting the cross correlation of the source function and its Hilbert 


transform) performs just as well as a 3-point match. The optimum offset parameter showed 


71 


to agree well with theory. Using the calculated amplitudes in SIFT instead of setting them 
to unity generally improves the results except when the amplitudes are totally in error. Thus, 
with amplitude variability being a major issue in real data, this approximation may be 
validated, in general. 

To be able to explain the good performance of the algorithm in the Barents Sea 
experiments, we note that the location of the experiment was more than 250 m deep, the 
signals were essentially low frequency, and the low frequency variant of the algorithms was 
used. The bottom was modeled as a hard bottom and attenuation was modeled as rough 
surface scattering, a reasonable approximation for this area. Both synthetic and real (SUS) 
sources were used. The source range was approximately 4.5 km, the source depth was 15 m 
and the receiver depth was 170 m for most of the analysis. Applying the results found 1n this 
thesis suggests that the larger depth was a crucial point in the performance of the algorithm 
in the Barents Sea analysis. In addition, shallow source depth and large grid spacing in range 
Cause any change in the shape of the matching point to be resolved in the same resolution 


cell. 
J. SIFT MATCHING RESULTS (HIGH FREQUENCY) 
Both the 1-point and 3-point high frequency SIFT matching algorithms were found 


not to be able to localize in the negative gradient environment beyond 4500 m. In addition, 


the 3-point algorithm was not be able to localize at ranges shorter than 2000 m. The 1-point 
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algorithm localizes correctly at ranges shorter than 2000 m but shows high sidelobes. 
Beyond 4500 m both algorithms show an area of high sidelobes around the actual source 
location. This is only seen for very high bandwidths (>250 Hz true bandwidth). 

The time domain arrival structure of the received signal for the 3-point match shows 
a large onset corresponding to the Hilbert transform component of the main arrival. This 
results in an offset in the autocorrelation which may partly explain the difficulties with the 
algorithm. Using ray code generated source signals shows good performance for both 
algorithms even at 250 Hz bandwidth. This suggests that it might even be a programming 
error since these results should experience the same difficulties. 

Both the 1-point and the 3-point algorithms show good performance for the iso-speed 
environment. Manipulating the individual amplitudes in the 3-point match does improve the 
results but only slightly. Other results, however. suggest that the add and compare approach, 
discussed earlier, may improve the 1-point and 3-point matching results at larger ranges. 
Some results suggest that both high frequency algorithms require even more bandwidth than 
the low frequency version of SIFT. Unfortunately, this might not be possible in practice. It 
is expected that there might be some better performance for the 3-point matching algorithm 
in the 2-4.5 km range. The more general n-point SIFT algorithm has not been used as it was 
expected to suffer even more from the limited bandwidth, and it was considered to be very 


difficult to tune the individual 6t's and amplitudes of the components correctly. 
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Figure 11 Low frequency SIFT results for the negative sound 
speed gradient environment using the 3-point SIFT algorithm 
and 250 Hz bandwidth. The scattered match results from working 
beyond the performance limitations. (a) Without add and 
compare approach applied. (b) With add and compare approach 
applied. 
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Figure 12 HF frequency SIFT localization results using a PE 
generated source, located at 5500 m range and 59.8 m depth and 
a receiver located at 40.2 m depth, for a negative soundspeed 
gradient environment. (a) 1-point SIFT algorithm, bandwidth 
500 Hz, no window applied. Add and compare approach enhances 
Pmeesesult. (b) 3-point SIFY algorithm, 500° HZ bandwidth wea 
Blackman window. Add and compare approach enhances the result. 
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V. SUGGESTIONS FOR FUTURE RESEARCH AND PRELIMINARY 
ANALYSIS OF OTHER MATCHING ALGORITHMS 


Directions for future research for single hydrophone Matched Field Processing is one 
of the goals of this research. Some of the suggestions for future research may be slight 
modifications and improvements to existing algorithms. Those suggestions and their 
motivation have already been mentioned with the description of the algorithms in Chapter 
3 and the analysis and results in Chapter 4. In this chapter, suggestions for more fundamental 


changes to existing algorithms and new algorithms are given. 


A. SOURCE DEPENDENCE 


The dependence on the source spectrum has been neglected throughout this research. 
If the source is a single, broadband, pulse-like transient, the specific details of the spectrum 
may not be critically 1mportant. A standard practice may be to filter out some portion of the 
signal spectrum over which all the processing will be done. An algorithm which 1s fairly 
insensitive to the spectral details 1s then needed. 

However, if the source is composed of two or more closely spaced pulse-like events, 
the separation of the various multipaths at the receiver location becomes a formidable task. 
At issue then is whether a localization algorithm can be effective when utilizing only 


snippets of the arrival (i.e., portions of the signal believed to contain multipaths from a 
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single pulse-like event). Presumably, this is true if enough single event multipath 
information is available for analysis. The degradation of the localization due to multipath 


information needs to be investigated further. 


B. MATCHING THE TIME-DOMAIN ARRIVAL STRUCTURES 


Earlier analysis which directly compared arrival structures indicated that this 
coherent approach is too sensitive to phase errors. The autocorrelation matching technique 
also suffers from phase errors in the phase differences between separate arrivals. This 
suggests that we first remove all phase information by taking the absolute value of the 
basebanded, time-domain arrival structure. In order to emphasize the smaller amplitude of 
later arrivals, it might also be beneficial to then apply a high-pass filter in the frequency 
domain to remove the mean and low-order trends from the time-domain, absolute value 
arrival structure. On this signal you can apply the standard autocorrelation matching or just 
apply a normalized cross-correlation with a replica to perform the matching. Careful design 
of the high-pass filter is required, however, due to the transient response of the filter in the 
time-domain (The large peak in the arrival structure due to the initial arrivals triggers a 
large transient response). Figs. 13 (a) and 13 (b) | show the absolute value of the arrival 


structure and.the demoduted signal. 
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A limited initial analysis has been done with some success. As the problem requires 
a specific filter which creates a ainiteam of phase and amplitude distortion a low 
coefficient Butterworth zero phase filter algorithm has been used. For the complete range 
in the iso-speed and negative gradient environment, a filter with 2 coefficients (effectively 
4 coefficients due to the zero phase filter algorithm) and a cutoff of 0.3 times the Nyquist 
frequency using a sampling frequency of 250 Hz gave reasonable results. (Nb. The cutoff 
is related to the bandwidth and the sampling frequency). As the parameters of the matching 
algorithm have not been tuned to optimum performance, results may change a little after 
tuning of the parameters (e.g. number of notched lag values). 

For autocorrelation matching, the baseline result using a negative gradient 
environment gives a footprint at 5500 m which 1s 80 cm in depth and approximately 40 m 
in range. Using a coherent match, the footprint becomes slightly larger in range and 
considerably larger in depth, approximately 120 cm in depth and 40 m range. Sidelobe 
levels are, however, considerably higher than for standard autocorrelation matching or 
coherent matching. Preliminary results using a ray code based template show that the direct 
match (normalized cross correlation) and the autocorrelation may perform very well for 
some source locations, while for other source locations performance may degrade severely 


due to high sidelobes. All footprints were found to be very small in depth. 
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It has been shown that the sample rate has a significant influence on the results. For 
the iso-speed environment, similar results have been found. The footprints for this 
environment are, however, larger. Combining this approach with a SIFT-like 1-point 
approximation algorithm did not perform well. This approach only relies on arrival times 


and completely neglects the phase and frequency for making the templates. 


w. PARAMETRIZATION OF THE TIME-DOMAIN ARRIVAL STRUCTURE 


Although originally not incorporated in the project, some analysis has been done with 
parametrization of the arrival structures and the replicas. As the arrival structures decay 
rapidly, it was suggested that they could be modeled as an auto regressive process. For the 
analysis, linear prediction coefficients were used and the Euclidian distance was used as a 
matching criterion. Results of the matching suggested that matching is possible without 
having to align the signals. Taking it a bit further would suggest that the same approach 
could be applied to one side of the autocorrelation function. 

A limited preliminary analysis of the results has been performed. The scale used to 
present the results is log(1/D) where D is the euclidian distance. The matching works like 
a signature. The baseline footprint is very small, 1 resolution cell (5 m x 40 cm). A limited 
analysis of the robustness has been done using a random uniform 20% disturbance of the 
amplitudes of the frequency bins of the PE results. The match will stretch in range to 2-5 


resolution cells and 2 resolution cells in depth. The algorithm works better for a low number 
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of coefficients. At a range of 5500 m for the negative gradient environment, performance 
is severely degraded using more than 10 coefficients which suggests that the description of 
the signals must not be too detailed. No phase perturbations have been tested. It 1s assumed 
that this might more severely degrade the performance. Using a PE generated source signal 
and ray code generated templates the algorithm still localizes but is severely degraded by 
high sidelobes. This preliminary analysis has only been performed in the 4500-6500 m 


range. 
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sidelobes. 
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D. AUTOCORRELATION MATCHING OF THE FREQUENCY DOMAIN 
RESPONSE 


The received signal shows a particular signature in the time domain as well as in the 
frequency domain due to the multipath interference. Note that the absolute value of the time- 
domain arrival structure described above is related to the autocorrelation of the frequency 


domain response through a Fourier transform by 


Gae= [ ROOR (Fe xax 
, (93) 
= | R(OR “(De -?™ dt. 


—OO 


A frequency-domain autocorrelation matching algorithm similar to the time-domain 
approach would then produce an ambiguity surface. The autocorrelation matching can now 


be formulated as 


A= | Gre()Gop Daf (94) 
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and we can again notch out values near zero “lag” to remove the mean and low-order trends 
of the time domain. Note that this technique is essentially equivalent to matching the 
absolute value of the time-domain arrival structures. 

Formulating the matching of the autocorrelation of the frequency response, however, 
has reintroduced the problem of defining an absolute reference time. In the definition of the 
replica autocorrelation, it is therefore necessary to introduce a free parameter to allow for 


searching in the time domain, 1.e., 
nll 2)* [Op “(t-t)e Mat 
“fo TRE )e Tae (95) 
=e Rf fPooe “(x fdf 
The autocorrelation matching would now look like 


Aor = i GG Gods: (96) 


The values used for the ambiguity surface are now the values of A(z for which A is 
maximum. Furthermore, this algorithm applies to both real and complex signals, an 


improvement over the original SIFT design. 
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It should be noted that, while we developed this idea independently, it was recently 
brought to our attention that other investigators have previously used an extremely similar 
approach with reasonable success (Nghiem-Phu et al., 1992). We believe they were the first 
to suggest using this approach for this problem. However, they did not apply any of the 
filtering techniques we have suggested and improvements to their original concepts may still 


be forthcoming. 
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VI. CONCLUSIONS 


Results of autocorrelation matching and approximate autocorrelation matching at the 
primary range of interest (5-10 km) proved to be highly inconsistent and sensitive to signal 
mismatch between ray-based methods and full-wave techniques. Due to our inability to 
generate data in the 2-4.5 km range for the negative sound speed gradient environment, only 
limited conclusions can be drawn from the results. The autocorrelation matching results 
appear to be bandwidth limited for longer ranges and limited by the sidelobe levels for 
shorter ranges. Results suggest that using the full 250 Hz bandwidth without windowing can 
lead to succesfull localization at ranges of more than 5 km. Footprint sizes are very small 
in depth (< 1 m) and moderately small in range (<<15 m). Taking into account phase 
differences in the autocorrelations improves the peak to sidelobe levels. Notching out values 
near zero lag improves the dynamic range and reduces the chance of invalid localizations 
(due to energy matching instead of shape matching) but reduces the footprint size. Results 
for the iso-speed environment are much better, as may be expected, and footprints are larger 
than for the negative gradient environment. 

The SIFT algorithm is very sensitive and more bandwidth limited than the 
autocorrelation matching. For the low frequency version of the algorithm, a first order 
description of performance limitations is suggested. The high frequency versions of the 
algorithm did not give consistent results for the negative gradient environment in either the 


750-1000 Hz or the 500-1000 Hz band but proved to work well in the iso-speed 
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environment. The high frequency versions seem to be even more bandwitdh limited than the 
low frequency versions and thereby severly reduces the range of applicability. Results 
suggest that localization can be performed at ranges less than 5 km for a bandwidth of 250 
Hz and a bottom depth of 100 m and moderate sound speed gradient environments. Using 
a large bandwidth ‘localization by sidelobes’ could be performed in some cases but with 
minimal resolution. Finally, additional improvements may be achieved for the negative 
gradient environment by taking into account the phase shift in the ‘half matching 
autocorrelation results, using no window (thereby using the full bandwidth) and extending 
the range to include the 2-4.5 km bracket. The same approach may improve the SIFT results 


where the multi-point approach may also be considered. 
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APPENDIX: SPATIAL BEAMFORMING AS A FORM OF MFP 


The autocorrelation matching in the frequency domain can be related to a generalized 
and normalized Bartlett beamformer. The main difference is the neglect of the off diagonal 
terms in the correlation matrix of the frequency domain received signal. This is closely 
related to the fact mentioned by Miller et al. (1996) that the correlation-based algorithms are 
likely to work better on noise-like signals than on sinusoidally signals. 

The approach presented here is an analogue of the spatial techniques known in 
beamforming and frequency bank techniques known in spectral estimation. In most 
applications of Matched Field Processing (MFP), an array of receiver hydrophones is used. 


The output of the beamformer for the standard Bartlett beamformer is then defined as 


B,=E{|D()P} =we R,We=w "E{p()p "()}w, (97) 


where R, is the matrix of the data p(t), and w is the weighting vector for a particular look 


direction. The data p(t) can be defined as 


P{t=a(O)s{H+n() — i=1,......,N, (98) 


where N is the number of elements in the array, 7 /t) is the noise, a;(@) is the array response 
for asource in the direction @, and s,(t) is the source function for a source in direction 6,. 
Usually this approach is based on the plane wave approximation, and the beamformer can 


be refined to multi-variance, multiple constraint or more exotic schemes. The step from 
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plane wave beamforming to MFP is made by including the knowledge of the environment 
and the propagation, usually by means of a propagation model, in the weighting vectors. 
This enables you to make estimates in more dimensions, usually a 2-D range/depth plane. 
Necessary for this generalized beamforming to work is that the description of the 
environment and the propagation is accurate enough, the environment has enough vaniation 
over the distances of interest, and the weight vectors at different points are different enough 
to form a consistent match. In MFP the weighting vectors are usually called replicas and the 
generalized beamformers are called processors. The most common are the Bartlett, the 


multi-variance and the multiple constraint processors, 
_.,H 
B,=w' Row, 
| Ap-l fl 
Byyy=\w Rew) (99) 


Bycy=e!|W4Rs'W| ‘ce, 


where W is the matrix of the replica vectors defining the multiple constraints around the look 


direction. The elements of c are defined by 


_.)H 
Cw <.. (100) 
where w, is the replica vector in the look direction and w,, are the vectors defining the 
constraints. 
The above spatial approach can be applied to a broadband time domain signal using 


one hydrophone instead of the spatial array. Two different approaches will be shown, one 
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in the time domain and one in the frequency domain. We define the preenvelope (or 


complex envelope of the signal) of the received pressure as 


P(t.) =as(t,)+n(t,) f=1.......N, (101) 


where @ 1s the complex response of a hydrophone to a source at location (r,, z, ) to a 


complex source at time t,. The Bartlett processor can be formulated as 


B,=w"R w (102) 


where w are the replicas. The problem with this time-domain approach is the lack of a time 
reference for the replicas. With an array, all the hydrophone locations are assumed known 
and provide a reference in space. The reference in time 1s not so easy to establish. Unless 
you can overcome this problem, this approach does not seem viable 

On the other hand, taking the second frequency-domain approach a proper absolute 


time reference can be established. Define the frequency domain received pressure as 


PL) =as(f)+n(7)). (103) 


Treating it the same way as the spatial array R, can then be defined by 


R= {P(f) B'(p} 


(104) 
=E{G S/S" Na" +a(f)n"(N}. 
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The generalized ‘beamformer’ or matching algorithms can be defined as 


B,=W"R Ww. (105) 


This approach does not encounter the reference problems an approach in the time domain 
would have. This approach could also be extended to other types of processors. The Bartlett 
processor is closely related to autocorrelation matching. If we define the autocorrelation 


matching in the frequency domain as 


Geyer) 


[SIG PD IGE 


A= (106) 


it is clear that this looks like a normalized Bartlett processor where the off diagonal terms 


have been neglected. 
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